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Abstract. In this work we concentrate on the evolution of the cluster population of the interacting galaxy M51 
(NGC 5194), namely the timescale of cluster disruption and possible variations in the cluster formation rate. We 
present a method to compare observed age vs. mass number density diagrams with predicted populations including 
various physical input parameters like the cluster initial mass function, cluster disruption, cluster formation rate 
and star bursts. If we assume that the cluster formation rate increases at the moments of the encounters with 
NGC 5195, we find an increase in the cluster formation rate of a factor of 3.01^2, combined with a disruption 
timescale which is slightly higher then when assuming a constant formation rate (£4 = 2.0 jli'f x 10 8 yr vs. 
l.Olg'ss x 10 8 yr). The measured cluster disruption time is a factor of 5 shorter than expected on theoretical 
grounds. This implies that the disk of M51 is not a preferred location for survival of young globular clusters, since 
even clusters with masses of the order of 1O 6 M0 will be destroyed within a few Gyr. 
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• 1— 1 1. Introduction 

X' 

$H ' The goal of this scries of papers is to understand the 
, properties of the entire star cluster population of the 
interacting spiral galaxy M51. These properties include 
the age and mass distribution of the cluster population. 
Additional properties are the survival rate of the clusters 
as well as any relations between the observed properties. 
These relations may be used to constrain cluster formation 
and destruction scenarios. 

In order to study the above properties, we exploit the 
large amount of HST broad-band archival data on M51, 
which covers roughly 50% of the observed surface area 
of M51, and covers a broad spectral range (UV to NIR). 
The large spatial coverage is necessary in order to obtain 
a large sample of clusters for carrying out a statistical 
analysis, and the broad spectral range allows accurate de- 
terminations of the individual cluster properties (Bik et 
al. 2003; Anders et al. I2004|l . A preliminary analysis of a 
subset of the M51 cluster population was carried out by 
Bik et al. (2003; hereafter Paper I) who introduced the 
method used to determine the cluster properties and de- 
rived the age and mass distributions of the cluster sample 
roughly 2 kpc to the North East of the nucleus. 
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Bastian et al. (2005; hereafter Paper II) extended the 
survey to include the entire inner ~ 5 kpc of M51, and 
found 1152 clusters, 305 of which had accurate size deter- 
minations. In that work we extended the age distribution 
analysis of Paper I and found evidence for a cluster for- 
mation rate increase ~ 50 — 70 Myr ago. This corresponds 
to the last close passage of NGC 5195 and M51 (Salo & 
Laurikaincn 2000). Additionally we found that 68 ±15% 
of the clusters forming in M51 will disrupt within the hrst 
~ 10 Myr after their formation, independent of their mass, 
so-called infant mortality. For the resolved cluster sample, 
we found that the size distribution (the number of clusters 
as a function of their effective radius) can be well fit by a 
power-law: N dr R oc r c ff _?J dr e g, with r\ = 2.2±0.2, which 
is very similar to that found for Galactic globular clusters. 
Finally, we did not find any relation between the age and 
mass, mass and size, nor distance from the galactic center 
and cluster size. 

In this study we focus on the evolution of the popula- 
tion of clusters in M51, in particular the timescale of clus- 
ter disruption and possible variations in the cluster for- 
mation rate. Cluster disruption of multi-aged populations, 
which excludes the galactic globular clusters, has been the 
subject of many earlier studies (e.g. Hodge l|1987|l for the 
SMC and Battinelli & Capuzzo-Dolcetta <|1991|l for the 
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Milky Way). In this study we will take into account mass 
dependent disruption, since the time needed to destroy 
half of the cluster population, which has been estimated 
in earlier work, will strongly depend on the mean mass 
of the sample and the lower mass limit of the sample. In 
addition, we here want to study the effect of variations in 
the formation rate, which is usually kept constant. 

Boutloukos & Lamers l|2003[) have developed a method 
to derive the disruption timescale based on the age and 
mass distributions of a magnitude limited cluster sample. 
They found that the disruption time of clusters in M51 
is a factor of 15 shorter than the one for open clusters 
in the solar neighborhood. Lamers, Gieles & Portegies 
Zwart (2005) showed that part of the difference can be 
explained by the difference in density of the cluster envi- 
ronment. They showed that the disruption time of clusters 
depends on the clusters initial mass and the galaxy density 
as tdis cx M° 62 p~®{ 5 , based on the results of iV-body sim- 
ulations. The disruption time of clusters in M51 was still 
about a factor of 10 lower than the predicted value. In this 
work we are particularly interested if a short disruption 
timescale can be mimicked by an increasing cluster forma- 
tion rate, and how the assumed disruption law influences 
the derived timescales. To this end we have generated ar- 
tificial cluster samples with parameterized global charac- 
teristics (e.g. time dependent cluster formation rates, dis- 
ruption laws, infant mortality rates, and mass functions). 
We then compare these models with the derived age and 
mass distributions of the cluster population of M51 to de- 
rive the best fit parameters for the population as a whole. 

The structure of the paper is as follows: In §[2]the ob- 
servations of the cluster population of M51 are presented. 
In §121 we investigate how the disruption time depends on 
different cluster and galaxy parameters. §01 describes the 
steps we will take in our models, where the details of the 
models we used to generate artificial cluster populations 
will be explained in § The results of the fits are given in 
§ El A discussion on the implication of the results is given 
in § The conclusions are presented in § |HI 

2. The observations 

2.1. Fitting the observed spectral energy distribution 

From archival HST broadband photometry we have de- 
rived the age, mass, and extinction of 1152 clusters in M51 
(Paper II) using the three-dimensional maximum likeli- 
hood fitting (3DEF) method. Details about the 3DEF 
method can be found in Paper I. In summary, the spec- 
tral energy distribution of each cluster is compared with 
cluster evolution models. In this case the GALEV sim- 
ple stellar population (SSP) models (Anders et al. [2003 
Schulz et al. I2002f) for solar metallicity and Salpeter IMF 
are used. For each age a series of different extinctions is 
then applied to the models and all combinations of age 
and extinction are compared to the data. The lowest x 2 
is kept and from the absolute magnitude at that age the 
mass is determined. Detailed tests of the accuracy and re- 



liability of the derived parameters are presented in Paper 
II. In the present work we further investigate the accuracy 
of our fitting method and we use the data set to develop 
a model that describes the global properties of the cluster 
system. 

2.2. The age and mass distribution of clusters 

The ages and present masses of the 1152 clusters are plot- 
ted in the top panel of Fig. ^ In order to be able to com- 
pare our observations with simulated cluster samples, we 
bin the data in logarithmic number density plots of the 
age vs. mass distribution. Clusters where counted in bins 
of 0.4 age dex by 0.4 mass dex (Fig.^ Middle). The result 
is illustrated in the bottom panel of Fig. ^ A few striking 
features can be learned from this diagram: 

1. a burst in the cluster formation rate (CFR) between 
50 and 70 Myr, corresponding with the most recent 
interaction with the companion galaxy NGC 5195; 

2. a short lived young population with ages < 10 Myr. In 
Paper II we found that ±68% of these young clusters 
will dissolve within 10 Myr, independent of their mass; 

3. evolutionary fading under the detection limit, which 
makes it harder to detect old low mass clusters. The 
increasing line in Fig. ^ shows how the 90% complete- 
ness limit in the F439W band (22.6 mag.) corresponds 
with different masses at different ages; 

4. an apparent increase in the mass of the most massive 
cluster with age. This is a binning effect: the older age 
bins span more time and therefore contain more clus- 
ters and the chance of finding a more massive clus- 
ter at older ages is higher due to the size of sam- 
ple effect (Hunter et al. 2003 s. We have tried to use 
the method of Hunter et al. ( 2003 ) to derive the clus- 
ter formation rate, but we found that the increase in 
the maximum mass is much to shallow. This proba- 
bly means that M51 has reached the maximum cluster 
mass (~ 10 6 M Q ) and then the relation of Hunter et al. 
( 2003J does not apply anymore, since the maximum 
cluster mass found at a certain age is then not deter- 
mined anymore by sampling statistics. This is the topic 
of a next study (Gieles et al. |2005 ) . 

We found that a large fraction (±68%) of the clus- 
ters younger than 10 7 yr dissolves independent of mass, 
probably due to the removal of the primordial gas (e.g. 
KroupaESlI Geyer & Burkert l2T)0TI Lada & Lada HSH - 
After these critical 10 7 years, the surviving clusters will 
dissolve due to the tidal field of the host galaxy and exter- 
nal perturbations from, for example, encounters with giant 
molecular clouds (GMCs). The disruption of clusters is in 
that sense a two step process. Here we will use the re- 
sulting age/mass distribution to study the disruption of 
clusters with ages larger than 10 7 yr., i.e. the second step 
in the disruption process. 
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Fig. 1. Ages and present masses of the 1152 clusters iden- 
tified in Paper II. Top: original data, where every point 
represents a cluster. Middle: same data as in top panel, 
with overplotted the grid used to bin the data. Above the 
dark line are the bins which are not affected by the detec- 
tion limit. Bottom: logarithmic density plot of the same 
sample, where dark regions represent more clusters. The 
right hand scale shows how the grey values correspond to 
the logarithm of number. The line in all three plots is the 
90% completeness limit (F439W = 22.6 mag). 



2.3. Artifacts introduced by the age fitting method 

We want to see whether the 3DEF method (§E3J, used to 
derive ages, masses and extinctions from the photometry 
introduces systematic artifacts. More important, could it 



affect our results of the disruption time or formation rate? 
For instance, are there systematically old clusters fitted 
with young ages or the other way around? 

The uncertainty in the derived ages, extinctions and 
masses from broad-band photometry is mainly caused by 
two effects: 

1. Differences between the real integrated colors of the 
clusters and the models, for example due to stochastic 
sampling of the stellar IMF which cannot be taken 
into account in the SSP models, or errors in the stellar 
isochrones used; 

2. Systematic errors introduced by the age fitting method 
and the applied selection effects 

Of course the first effect cannot be corrected for, un- 
less we are able to compare photometric determined ages 
directly with spectroscopically determined ages, which is 
unfortunately not feasible for a whole population of clus- 
ters. In addition, although spectroscopically derived ages 
are more accurate, they are hampered by their own prob- 
lems (Brodie et al. 11998(1 . In addition, the age derivation 
will dependent on the choice of the adopted SSP mod- 
els. Variations in the metallicity and the IMF will af- 
fect the derived ages. A detailed comparison of the data 
with models of different metallicity is carried out in Paper 
II. An earlier study of the age distribution of M51 used 
Starburst99 models (Bastian & Lamers l2003[> . and they 
found a very similar age distribution. 

The second effect can be quantified with the use of ar- 
tificial cluster populations. Earlier studies (e.g. Anders et 
al. 120041 de Grijs et al. l2()05|l have already shown the im- 
portance of using a long wavelength baseline (U to NIR) 
for age dating young clusters. Here we will make an at- 
tempt to quantify possible systematic errors introduced 
by the age-fitting method and see whether we can correct 
for them or not. 

To quantify the artifacts introduced by the fitting rou- 
tine, an artificial cluster sample including simulated ob- 
servational errors and extinction values is generated and 
fitted with the same fitting procedure as used for the data 
(§H3J. We start with a sample of clusters equally spread in 
log(Age/yr) and log(M/M Q ) space. In total 201 time steps 
between log(Age/yr) = 6 and log(Age/yr) = 10 and 161 
mass-steps between log(M/M Q ) = 2 and log(M/M ) = 
7 were generated. The GALEV models have log(Age/yr) 
= 6.6 as youngest model, so clusters with younger ages 
were given that age. The magnitudes as a function of 
age and mass where taken from the GALEV SSP models. 
Observational uncertainties were applied as a function of 
magnitude as was done in Paper II: the observed errors 
in the magnitudes of clusters in M51 can be well approxi- 
mated by Amag A = io d i+ d 2xmag A Thc values for ^ and 

di for the different filters are results from analytical func- 
tions fitted to the observed errors and magnitudes and are 
given in Table 4 of Paper II. Ideally, we then apply the 
same extinction to the model clusters as the M51 clus- 
ters have. Unfortunately the only information we have is 
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the extinction we have measured, which of course could 
already be polluted with artifacts. To get an estimate of 
the uncertainty in the measured extinction, we start with 
a sample of clusters with no extinction applied. When we 
fit this population with the 3DEF method, we find that 
20% of the sources is fitted with some extinction. This is 
quite a large number, but fortunately 90% of these sources 
have extinction values lower than E(B — V) = 0.1 mag. 
The maximum E(B — V) found is 1 mag. 

The next step is to apply an extinction model close to 
what we observe. To this end E(B — V) extinction val- 
ues were chosen randomly from a Gaussian distribution 
centered at with a = 0.10 for clusters younger than 
log(Age/yr) = 7.3 and a = 0.05 for clusters older than 
log(Age/yr) = 7.3. The values for a agree with the value 
we found for the mean extinction in Paper II. There we 
found that these values are the average extinction for these 
two age groups. The higher extinction for young ages is 
caused by the presence of the left-over dust around the 
cluster. Negative extinctions were set to 0, resembling the 
extinction distribution of the data where half of the clus- 
ters had E(B - V) = (See Fig. 8 of Paper II). An age 
dependent maximum extinction was applied of the form: 
E(B - V0max(*) = 5 - 0.5 * log(i). This is a little bit lower 
than the observed maximum extinction, but we know that 
some of the observed high values could be caused by wrong 
fits. This still resembles the observed extinction behavior 
quite well. The resulting magnitudes are than cut off at 
our completeness limits in each filter. In this way we cre- 
ated the spectral energy distributions of a large artificial 
cluster sample with age, mass and extinction known for 
each cluster. These were fitted with the 3DEF method 

(§EU). 

The result of the fitted simulation is shown in Fig. [21 
A direct comparison with the observed age-mass diagram 
of M51 clusters (Fig. [I] Top) shows that there are features 
present in the data which are not visible in the fitted sim- 
ulation. For example, there is a gap at 6.9 < log(Age/yr) 
< 7.1 in the M51 cluster sample which seems to appear 
at slightly higher ages in the fitted simulations (7.1 < 
log(Age/yr) < 7.2 ). This suggests that the artifacts in 
the data are not only caused by our applied selection ef- 
fects or our age-fitting technique. In the top left panel 
of Fig. we show the fitted age versus the input age for 
the simulated cluster sample. A large number of clusters 
with wrong ages are fitted with an log(Age/yr) ~ 7. We 
found for 87% of the modeled clusters that the fitted age 
was the same as the input age within 0.4 dex (Fig.|21 Top 
right). For the mass 97% was fitted correctly within 0.4 
dex (Fig. Bottom right) and 92% of the extinction val- 
ues where fitted back within 0.05 mag. (Fig. Bottom 
left). We have to realize that the strength of the artifacts 
depend on the number of input clusters at each age and 
mass bin. We have not attempted to match the observa- 
tions in this stage, since we are only interested in relative 
errors. For example, the horizontal spur at log(Age/yr) ~ 
7 in Fig. |3] (Top, left) is populated with clusters with in- 
put ages up till a few times 10 9 yr. The number of clusters 




Fig. 2. Age/mass diagram of artificial sample after fit- 
ting with the method described in § 12.11 The deviations 
from equally spaced dots in log(Age/yr)/log(A^f/M Q ) are 
caused by the applied observational errors and the fitting 
routine. 
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Fig. 3. Result of fitting an artificial cluster sample with 
the 3DEF method. Top left: The fitted age is shown ver- 
sus the input age for each cluster. Deviations from the one- 
to-one relation are caused by photometric errors which are 
applied to the input sample and misfitted extinction. Top 
right: The percentage of clusters fitted with the same age 
as the input value plus some deviations, as a function of 
this deviation. Bottom left: similar but then for extinc- 
tion. Bottom right: similar but then for the mass. 



with that age in our M51 sample is very low (See Fig.^ 
Top). 
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2.4. Correcting the data for fitting artifacts 

We will try to correct the data for possible artifacts, by us- 
ing the (systematic) deviations found in § 12.31 Correcting 
the observed ages based on the absolute numbers deviating 
from the one-to-one relation is not useful, since the num- 
ber of clusters that were used as input at each age and 
mass differs from the observed number. From the input 
sample we can derive how many clusters are (systemat- 
ically) fitted with wrong ages and masses. Let the total 
number of bins in age/mass space be K. Here we define 
the number of bins as the number of bins which are not 
affected by the detection limit (See Fig.f^ Middle panel). 
The number of clusters found in each bin is the sum of the 
contribution of clusters from all bins to this one, where the 
majority will be from the bin with the same input age and 
mass. When we write the number of clusters in each bin as 
a vector with if -entries, the fitted number of clusters can 
be written as a matrix multiplication of all contributions 
times the input number of clusters 





( 








obs V 



^11 

C21 



Cl2 
C22 



c 



K2 



c 



KK 



\ 








1 


\n k ) 



(1) 



where N^.j is the number of clusters generated in bin 
j, N b s ^ is the number of clusters fitted in bin i and Cy 
is the contribution to bin i from bin j. Tests have shown 
that the best results are acquired when only taking bins 
into account which are not affected by the detection limit 
(see Fig.f^ Middle panel). From the simulated and fitted 
sample we can derive the values for Cy for all combina- 
tions of i and j. All values on the diagonal of C (i.e. CV, 
where i = j) are close to 1 since most clusters are fitted 
with the same age, extinction and mass. All other values 
are or between and 1. When we know the matrix C, 
the inverse can be used to correct the observations for 
systematic 3DEF fitting artifacts. 



iV intr = C- 1 x N 



bs 



(2) 



where Ni ntT is the vector with the intrinsic number of 
clusters, C" 1 is the inverse of the contribution matrix as 
defined in Eq. ^ and 7V obs is the vector with observed 
clusters. When we correct the observed vector and plot 
it in a 2D age-mass diagram again, divide the corrected 
and uncorrected observation, we can see where deviations 
take place. Fig. 0] shows the ratio of corrected over un- 
corrected observations. The number of clusters in the age 
bins at log(Age/yr) = 7.2 and 8.0 has been lowered with 
about 15%. The observations, however, show a gap at 
log(Age/yr) = 7.2 (Fig. [IJ Top). So, the underestimation 
of clusters in that age bin is not caused by our fitting 
routine. 

The corrected observations based on Eq.|2are shown in 
Fig. El The burst at between 50-70 Myr is less pronounced, 
but still present. The differences with the uncorrected ob- 
servations (Fig. 2] Bottom) are small, and we therefore 
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Fig. 4. Ratio of the corrected observations over the uncor- 
rected observations. The corrected data was calculated us- 
ing Eq.QJand Eq. [3 Light regions indicate where less clus- 
ters are found back by the fitting procedure. Dark regions 
indicate where more clusters are fitted than inputted. 
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Fig. 5. Corrected age and mass distribution. The raw data 
from Fig. ^ (Bottom) has been multiplied with the inverse 
of the contribution matrix (Eq. ^) 



conclude that our age-fitting method (3DEF) and our ap- 
plied selection affect is not severely affecting our age-mass 
diagrams in a systematic way. Especially, there is no large 
systematic shift from old to young clusters or the other 
way around. We therefore conclude that we can use the 
uncorrected data as well as the corrected data to compare 
with the synthetic cluster populations in § El In § 16.21 we 
will show that both the corrected as the uncorrected data 
give the same results when fitting the analytical models 
to the data. 



3. Exploration of the parameters which determine 
the disruption time 

3.1. Cluster initial mass 

If cluster relaxation drives the evaporation of clusters, 
then the more massive clusters live longer than their low 
mass counterparts. Boutloukos & Lamcrs ( 2003) have pro- 
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posed an empirical way to determine the dependence of 
the cluster disruption time on the initial cluster mass, as- 
suming a power-law dependence of the disruption on the 
cluster mass 

i dis = t 4 (M t / 10 4 M o )t (3) 

where t± is the disruption time of a 10 4 Mo cluster, Mi is 
the initial mass of the cluster and 7 is a dimensionless in- 
dex. The value of 7 was determined by measuring the slope 
of the age and mass distributions of a cluster population. 
Their mean value for 7 based on four different galaxies was 
< 7 >= 0.62±0.06. In a recent study by Lamers, Gieles 
& Portegies Zwart (J2005:) these observational results are 
compared with results of A^-body simulations. The value 
of 7 can be explained by tidally driven relaxation and was 
confirmed to be 0.62 by TV-body simulations. The expla- 
nation for this is that the disruption time of a cluster in 
a tidal field, depends on the relaxation time (t re i) and the 
crossing time (t cr ) of the cluster as (Baumgardt 2001): 

t dis oc t* el x tl~ x (4) 

Using the expression for t Te \ and t cr from Spitzer (^987), 
this implies: tdis oc (3{N/ In N) x . Baumgardt & Makino 
(2003) found two combinations of (3 and x, depending 
on the concentration of the clusters. Lamers, Gieles & 
Portegies Zwart ( 2005 ) showed that for both combinations 
tdis could be well approximated with tdis oc TV 62 oc M 62 , 
where N is the number of stars in the clusters and M is 
the total mass of the cluster. King ( 1958 ) already had the- 
oretical arguments to expect that the lifetime of clusters 
should depend on the mass as tdis M 2 ' 3 . The agreement 
between observations and TV-body simulations was first 
noted by Gieles et al. ( 2004). The physical background of 
why the disruption time does not scale directly with the 
relaxation time is given by Fukushigc & Heggic ( 2000 ) . 

3.2. Cluster radius 

Young clusters are not only affected by the external tidal 
field of the host galaxy, but they also undergo shocks from 
(giant) molecular clouds. Both these processes shorten the 
lifetime of clusters. For both cases the radius of the cluster 
is an important parameter in determining how fast the 
cluster will disrupt. However, both processes depend very 
different on the radius. From Eq.0]and the expression for 
the relaxation time and crossing time it follows that for the 
tidally driven relaxation the disruption time depends on 
the radius as tdis oc r 3 / 2 . Larger clusters live longer since 
they have a longer relaxation time, so it takes more time 
for stars to reach the tidal radius and leave the cluster. 
Spitzer ( 1958) has shown that the time needed for a cluster 
to get unbound due to external shocks, relates to the half 
mass radius of the cluster as: t s h oc r^ 3 , so here larger 
clusters live shorter (for isolated clusters). 

To see whether the radius of a cluster is an important 
parameter in disruption, we use the radii measurements 
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Fig. 6. Number density plot of r e g vs. age. The decreasing 
number going towards larger radii shows the power law 
behavior of the radius distribution. The right hand side 
shows how the grey values correspond to the logarithm of 
number of clusters. 



of Paper II. There we measured the projected half light 
radius (or effective radius) r e ff, which relates to the half 
mass radius as: r cff = 3/4 r h (Spitzer 1987). We make a 
number density plot of r e g vs. age for all clusters (Fig.EJ). 
There are clearly no old clusters with large radius, while 
the opposite is expected due to the size-of-sample effect 
(Hunter et al. l2003f) . This suggests that large clusters are 
disrupted preferentially. This suggests that shocks may be 
the dominating disruption effect. However, when a large 
fraction of the clusters is removed, independent of radius, 
the upper radius would also go down. This is a result from 
number statistics: less clusters in a power law distribution 
will result in a lower maximum value. So what really mat- 
ters here is whether the slope of the radius distribution 
changes in time or not. In Paper II it is shown that the 
cluster radius distribution of M51 is N(r) dr oc r~ v dr, 
with ?y = 2.2 ± 0.2. To see how the slope of the distri- 
bution depends on age, we divide our cluster sample in 
young (log(Age/yr) < 7.5) and old (log(Age/yr) > 7.5). 
Dividing the sample at log(Age/yr) = 7.5 yields two sam- 
ples of more or less equal size, which gives similar errros in 
the fit to both distributions. When we determine this in- 
dex rj for only young clusters we find 77 = 2.0 ± 0.4 and for 
old clusters we find 77 = 2.5 ± 0.6, which is very similar for 
the value found for the globular clusters in our Milky Way 
(77 = 2.4 ±0.2, Paper II). Although the radius distribution 
seems to get steeper with age, the errors are too large to 
place a strong constraint on this. We therefore do not take 
the radius into account as a free parameter when mod- 
eling the cluster disruption. Futher studies of M51 with 
higher resolution, for example with the Advanced Camera 
for Surveys (ACS), could shed light on how the radius of 
clusters affects the lifetime. 
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Fig. 7. Ratio of the number of clusters (N) between 3-5 
kpc and 1-3 kpc as a function of age. Overplotted is a 
model predicting this ratio for two disruption times dif- 
fering a factor of 1.8, based on Eq|SJ 

3.3. Distance to the galactic center 

Lamers, Gieles & Portegies Zwart (2005) and Baumgardt 
& Makino |2003§ have shown that the disruption time is 
expected to depend on the galactocentric distance of the 
cluster, the orbital velocity in the galaxy and the ambient 
density of the galaxy as 

idis « R G /V (5) 

<* Pa^b (6) 

where Rq is the distance to the galactic center, V is the ro- 
tational velocity of the cluster in the host galaxy at that 
distance and p a mb is the ambient density of the galaxy 
at the location of the cluster. The relation with the am- 
bient density holds only when a logarithmic potential is 
assumed. Since we are dealing with a disk galaxy, it is not 
so straightforward to derive the ambient density from the 
spherically symmetric logarithmic potential. Therefore, we 
prefer the relation with the galactocentric distance and the 
velocity (Eq.EJ). With Eq. 0we are able to estimate if we 
would be able to observe a difference in disruption time at 
different locations in the galaxy. When we look at clusters 
between 1 and 3 kpc and at clusters between 3 and 5 kpc, 
the average value of Rq goes up with a factor of 2. The 
rotational velocity of M51 increases from 200 kms -1 to 
225 kms^ 1 (Rand lT^Mj) . So from Eq.Elwe expect the dis- 
ruption time in the two samples to be different by a factor 
of 1.8. In Fig. [7| we plot the ratio of the number of clusters 
at different ages for the outer region (3-5 kpc) and inner 
region (1-3 kpc). Overplotted is the predicted ratio using 
Eq. |H1 for disruption timescale that is different by a factor 
of 1.8. Just as for the radius dependence we see that the 
distance to the galactic center plays a role, but the data 
is not sufficient to include it in our analysis. 

In conclusion, we see evidence for radius and galac- 
tocentric distance dependent disruption, but the noise is 
too large to include these parameters in the models. The 
mass of the cluster is the most dominating parameter in 



the determination of disruption time and in the remainder 
of the study we will only use the mass dependence as a 
parameter we will vary in the models. 

4. Input parameters for modeling the cluster 
population of M51 

So far, analytical models for finding the cluster disruption 
time have assumed that clusters were formed with a con- 
stant CFR, as is probably the case for Galactic open clus- 
ters (Boutloukos & Lamers 120031 Battinelli & Capuzzo- 
Dolcetta I1991JI . Lamers et al. ()2005|l predicted the age 
distribution of open clusters. In the case of M51, we have 
age and mass information available for each cluster, so 
predictions can be done for age and mass. In addition, as- 
suming a constant CFR for M51 might be an oversimpli- 
fication of the situation, since the galaxy is in interaction 
with NGC 5195. In the next sections we explore a broader 
parameter space. 

Since there are strong arguments to believe that the 
mass dependence of the cluster disruption (7 = 0.62) is 
constant (Lamers, Gieles & Portegies Zwart 2005;. we 
start by varying only the constant £4, to be able to com- 
pare our results with clusters gradually loosing mass with 
the instantaneous disruption assumption (Eq. [3J| results 
of Boutloukos & Lamers (2003). Next, a two dimensional 
parameter search for 7 and £4 is performed, to verify the 
assumed value for 7 and to study the dependence of t\ on 
the value of 7. 

When we have a first estimate of the disruption time, 
we will study how this value changes when we assume that 
the CFR has been increasing during the last Gyr or con- 
tains bursts at the moments of encounter with NGC 5195. 

5. Analytical model for generating a cluster 
population 

5.1. Setting up a synthetic cluster population 

The synthetic cluster populations will be created in a sim- 
ilar way as in § 12.31 This time however we want to include 
realistic input physics, like the cluster IMF and differ- 
ent formation rates, so creating clusters equally spaced in 
log(Age/yr) and log(M/M0) will not be adequate. When 
creating clusters with a realistic CIMF, the number of 
clusters needed to fully sample the CIMF up to 10 Gyr 
ago is too high. Therefore each cluster was assigned a 
weight depending on the initial cluster mass and its age 
(w(t, Mi)), proportional to the expected number of clus- 
ters formed at each age and mass. The weight is a function 
of age and mass, scaled such that the youngest most mas- 
sive cluster has a weight of 1 

w(t, Mi) = (t/t min ) x (a - 1) x {Mi/M^) 1 -" (7) 

where w(t, Mi) is the weight assigned to a cluster with 
age t and mass Mi , M max is the mass of the most massive 
cluster in the simulation, t mul is the age of the youngest 
cluster in the simulation and a is the slope of the mass 
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function. When a is chosen 2, i.e. N{M) ~ A/~ 2 , the 
weight depends on age and mass simply as: w(t,Mi) oc 
t/Mi. When the simulated clusters are binned, the weights 
of the clusters are counted, yielding a realistic log(Age/yr) 
vs. log(M/M Q ) diagram similar to Fig. ^ (Bottom). The 
advantage of using points spread equally in log(Age/yr) 
and log(M/M©) with weights assigned, is that the number 
of points per bin is constant and that it is very easy to 
create a lot of populations with different formation rates, 
disruption timescales etc. in a short time. 

In our case the clusters are given a weight such that 
after binning the CIMF has a slope of a = 2.1 as found for 
M51 (Paper I) and the Galactic open clusters (Battinelli 
et al. I1994f) and using different formation and disruption 
scenarios f§ 16.21 5 IrTljl . 

5.2. Including stellar evolution and cluster disruption 

Baumgardt & Makino (2003) have shown that stellar evo- 
lution (SEV) is an important contributor to the disso- 
lution of young clusters, especially for clusters with low 
concentrations. They also confirmed that clusters dis- 
solve with a power-law dependence of their initial mass 
as tdis ~ M? i where 7 = 0.62, in agreement with the em- 
pirical determination by Boutloukos & Lamers ( 2003 ) . In 
the latter study instantaneous disruption after the disrup- 
tion time was assumed as a first approximation and they 
found that the typical disruption time (£4, see Eq.[3J) varies 
a lot for different galaxies. In a recent study (Lamers et 
al. |7005) it was shown that there is a simple analytical 
description of the mass of a cluster as a function of time. 
It takes into account the effect of mass loss due to stellar 
evolution, based on the mass loss predicted by the GALEV 
SSP models (Anders et al. I2UU51 Schulz et al. EUOljl and 
cluster mass loss due to the tidal fields. The mass of the 
cluster as a function of time can be well approximated by 

M p (*) = ((M iMse v) 7 -7f) 1/7 (8) 
to 

where M p (t) is the present mass of the cluster as a func- 
tion of its age, is the initial mass of the cluster, 
^ sov = M p (t)/Mi is the fraction of remaining mass af- 
ter mass loss due to stellar evolution has been subtracted 
and to relates to t± as t± = to x IO 47 . The mass as a func- 
tion of time, according to the analytical formula, agrees 
perfectly with the predictions following from A-body sim- 
ulations. Lamers et al. (J2D.05 ) have also shown that with 
this analytical model the age distribution of galactic open 
clusters can be explained very well. 

6. Fitting observed age-mass distribution to 
predictions 

6.1. Determining reduced x 2 values from 2D fits 

Artificial cluster samples with realistic input physics (e.g. 
a CIMF, cluster disruption, bursts etc.) can now be gener- 



ated and compared with the observed age and mass num- 
ber density distribution. 

After calculating the analytically generated cluster 
population, the model is binned into number density plots 
in the same way as the observed data (see § I2.2fl taking 
into account the weights. In order to compare the simu- 
lated (2D) age-mass density plots with the observations, 
we use the Poisson Probability Law (PPL) introduced by 
Dolphin & Kennicutt (2002) for similar purposes 

iv 

PPL = 2 V rrii - n t + n t In — (9) 
^ rrii 
i=0 

where N is the number of bins, rrii is the predicted number 
by the analytical model in bin i and rii is the observed 
number of clusters in bin i. The value of PPL is similar 
to the x 2 , m the sense that lower values imply better fits. 
We will always divide the PPL value by the number of 
bins minus the number of degrees of freedom, which is 
equivalent to the reduced x 2 , so the xt- We will refer to 
xl when we discuss results of fits. 

6.2. Determining the cluster disruption time assuming 
a constant formation rate of clusters 

To determine the typical cluster disruption time, £4, de- 
fined in § [Sj we generate a cluster sample with a constant 
CFR and then calculate the cluster masses as a function 
of age according to Eq. [S] for various values of i 4 . Here 
we are interested in the disruption time of clusters that 
have survived that first 10 7 yr in which the natal cloud is 
being removed by stellar winds, therefore we exclude the 
youngest age bin in the fits. Fig. |H1 shows a clear xt min- 
imum around <4 = LO^oJ x 10 8 vr > where the upper and 
lower errors are defined by xl, accept ^ X^min + 1 > which is 
equivalent to the 1 a error. In addition, we have fitted the 
same models but then corrected for age-fitting artefacts 
(§ 12. 4|) . The shape of this xt curve is the same as for the 
raw data, though the values are higher. This shows that 
the uncertainties of our age-fitting method do not alter 
the value found for the disruption timescale. 

To see how the value of £4 depends on the value of 7, 
we simulate a grid of cluster populations and vary £4 and 
7. A 2D xt plot is shown in Fig. The minimum is at 
7 = 0.65±g;Jg and U = l-0±S:fi x 108 y r > agreeing very 
well with the value of 7 = 0.62, which was stated earlier 
based on theoretical arguments and other observational 
results. The plot also shows that there is a diagonal bar- 
shaped minimum for different combinations of t$ and 7. 
One could argue that there could be multiple combina- 
tions possible, which will yield a somewhat higher value 
for £4. The fit however is very sensitive for the choice of 
bin size when varying two variables. We excluded the mass 
bins higher than 5 x 10 5 M Q , since we probably deal with a 
truncation of the mass function. If sampling effects would 
determine the upper mass at different ages (Hunter et 
al. I2003|l . the maximum mass should increase much more 
than we observe in the top panel of Fig. ^ This effect 
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Fig. 8. xt values for different disruption times. A clear 
minimum is visible at t 4 = 1.0±g;| x 10 s . The dotted line 
indicates xt m i n + 1- The dashed line shows a fit to the 
data after correcting for age-fitting artifacts (§ I2.4fl . In 
this simulations we have chosen 7 to be 0.62, based on 
theoretical arguments (§EJl- 
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Fig. 10. Illustration of the applied CFR increase in § 16.31 
for two different models. Model 1.: The CFR is taken to 
be constant before 1 Gyr ago and then increases linear in 
time until t = 0. Model 2.: The CFR increases with a step 
at the two moments of encounter with NGC 5195. Here 
the height of the step is the variable. 




14 [yr] 



Fig. 9. Two dimensional xt plot of 7 vs. £4. The cross and 
the shaded area indicate the region where xt m i n < xt < 
xl,min + 1 in *4 and 7. 



6.3. The effect of an increasing cluster formation rate 

Since NGC 5195 is probably bound to M51 and there- 
fore slowly falling in (Salo & Laurikainen 2000), one could 
argue that the short disruption timescale found in § 16.21 
is actually caused by an increasing cluster formation rate 
(CFR). Bergvall et al. (2003) have shown that interact- 
ing galaxies such as M51 (i.e. non-merging), can have an 
increased star formation rate of the order of a factor of 
2-3. We therefore model different cluster populations with 
increasing CFR(t) rates of various strengths, where we as- 
sume that an increasing star formation rate results in an 
equally large increase in the CFR(i). We study two dif- 
ferent models with increasing CFR: 1.) a linear increasing 
CFR starting 1 Gyr ago (§ I6.3.1|I : 2.) a CFR that increases 
with bursts at the moments of encounter with NGC 5195 
(S 16.3.2(1 . Fig. 1101 gives a schematic illustration of how the 
CFR varies with time for the two models. 



makes the mass function steeper above log(A//M Q ) ~ 5.3 
and therefor that region in the age/mass diagram is not 
suitable to fit the (sensitive) mass dependent disruption. 
An alternative way to measure 7 would be to measure 
the slopes of the age and mass distribution separately, as 
was done in Boutloukos & Lamcrs ( 2003 ) . We fitted these 
slopes and found the same value for 7 as for the 2D fit 
shown in Fig. El Again, for the mass, we do not include 
the high mass end for similar reasons as mentioned be- 
fore. This method is less sensitive for the choice of bin 
size, since we can fit the slope of the age and mass dis- 
tribution independent of the value of the disruption time. 
We choose to include the result of the simultaneous fit 
of ti and 7, because it illustrated nicely how these two 
variables relate. 



6.3.1. Linearly increasing cluster formation rate 

In the linear model the CFR(t) starts to increase 1 
Gyr ago, which is before the moment of the early close 
encounter with N GC 5195 (400-500 Myr ago, Salo & 
Laurikainen 2000). We expect the CFR to start increas- 
ing before the moment of the closest encounter, since the 
two galaxies are already interacting before the first peri- 
galactic passage. We calculate models with different CFR 
increases and disruption times. We plot the xt values for 
various values of CFR(i = 0)/CFR(t = 10 9 yr) and t 4 in 
top panel of Fig. ITT1 

The minimum xt value is at £4 = 2. Q*i\ x 10 8 yr and 
an increase m the CFR of 7.0±^ x . For small values of 
ti the equal xt lines are vertical. This can be explained 
by the fact that if the disruption time is short, no fin- 
gerprints of the ancient formation rate are present in the 
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current population. They are simply erased by disruption, 
the reason that the equal xt contours are circular around 
the minimum, is that the disruption of clusters depends 
on the mass of the clusters (Eq. [3J|, unlike an increase of 
the formation rate. 



6.3.2. Cluster formation rate with bursts 

An alternative formation scenario would be that the CFR 
increases with a burst at the moments of encounter with 
NGC 5195 and then an exponential decay in the CFR (see 
model 2 in Fig. ll(J|) . We choose the moments of increase at 
t = 7 x 10 7 yr and t = 5 x 10 8 yr ago, based on the results 
of Salo & Laurikaincn ( 2000 ) and the typical decay time 
of the burst is 10 s yr (Paper II). The CFR step and t± are 
varied in different models. The bottom panel of Fig. ITD 



+2.3 
-1.1 



X 10 s 



shows that the lowest xt value is at £4 = 2.0 
yr. This is a factor of 2 higher than when the increasing 
CFR is not taken into account, but it is the same value as 
was found for the linear increase in the CFR. The value 
is still a factor of 5 lower then predicted by ./V-body sim- 
ulations (Baumgardt & Makino 120031 Lamers, Gieles & 
Portegies Zwart I2()05fl . The best value for the increase in 
CFR at the moment of encounter is 3.0l^2- The latter 
value agrees very well with what is generally observed for 
the increase in star formation rate of interacting galaxies 
(Bergvall et al. 2003). Since one of the bursts is clearly 
observed and a linearly increasing CFR is not so physi- 
cal, we prefer Model 2 above Model 1 . In the next section 
we will compare several properties of this model with the 
observations. 



6.4. Comparision between the best fit model and the 
observations 

We show a direct comparison between the age-mass dia- 
grams of the best fit model fS l(5.3.2|l and the observations 
in Fig. The densities are scaled such that the total 
number of simulated clusters equals the total number of 
observed clusters (1152). A few bins in the observations 
are empty and not empty in the simulations. The reason 
for this is that the simulated cluster sample containts bins 
with values smaller than 1. Apart from this, the general 
trend of grey values in this 2D plot is very similar in both 
cases. 

Another interesting property of the observations is the 
formation rate. In Paper II we showed the number of clus- 
ters at different ages for different mass cut offs. For clus- 
ters with masses higher than 10 4 7 Mq we get a realistic 
impression of the cluster formation rate. This is because 
we are complete until 1 Gyr for these masses (see top 
panel of Fig.^) and because the most massive clusters are 
not affected by disruption that much. In Fig. 1131 (Toy) we 
show the number of clusters in different age bins for the 
observations and the best fit model. The general trend of 
the observations is followed very well by the model. A bet- 
ter way to show the formation rate is to divide each age 




M [yr] 




M [yr] 

Fig. 11. xt values for different combinations of CFR in- 
crease and disruption time (£4). The top panel shows the 
results of Model 1 where the CFR increases linearly be- 
tween t = 1 Gyr and t = present time. The bottom panel 
shows the result for Model 2 where the CFR increases 
with steps at the moment of encounter with NGC 5195. 
The cross and the shaded area indicate the region where 

Xi/,min Xu ^ Xi/.Tuin ^ 



bin by the width of the bin. Then we get the number of 
clusters formed per unit of time (Myr). This is shown in 
the bottom panel of Fig.^J In this figure the over-density 
of young clusters (log(Age/yr) < 7) is more obvious and 
the burst at 7 x 10 7 year is better visible. The first burst 
of cluster formation (5 x 10 s years ago) is not visible any- 
more, since clusters with these ages are already affected 
by the (short) disruption time. This reinforces that it is 
very hard to detect variations in the cluster formation rate 
when the disruption time is that short. The largest differ- 
ence is seen for the bin with log(Age/yr) = 7 and 8.25. 
The model predicts in these bins more clusters than are 
observed. This can be explained by fitting artefacts which 
yield an (unphysical) underdensity of clusters f§ l2.3[) . The 
model is still within the 3 a error of the observations how- 
ever. 
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Fig. 12. Comparison between the observed (Top) and the 
modeled (Bottom) age vs. mass number density plots. In 
both plots the 90% completeness limit of the F439W band 
is indicated with a line. The right hand side shows how 
the different grey values correspond to the logarithm of 
number. The total number in the simulations is scaled 
to the total number of observed clusters above the 90% 
completeness limit (1152). 

7. Implication of the derived disruption time 

We have shown that the clusters disruption time for a 
typical cluster with mass of 10 4 Mq is around 10 8 years in 
M51. When increasing formation rates are taken into ac- 
count the disruption time increases with a factor of 2. This 
is significantly longer than Boutloukos & Lamers (2003) 
found for clusters in a smaller region of M51 (£4 = 4 x 10 7 
yr). This could be because they did not seperate the dis- 
solution due to infant mortality rate from the evapora- 
tion by the tidal field from the galaxy. Clusters with ages 
younger than 10 7 year are not taken into account in this 
study, since they are affected by the dissolution due to the 
removal of primordial gas. Theoretical predictions show 
that clusters of 10 4 M Q in a tidal field of the strength of 
M51 should have an disruption time of about 10 9 years 
(Lamers, Gieles & Portegies Zwart 120051 Baumgardt & 
Makino 2003) . This value is found from observations of 
clusters in the solar neighborhood (Boutloukus & Lamers 
2003). What causes the clusters in M51 to dissolve about 
5 times faster than predicted? A few effects that have not 



og(Age/yr) 

Fig. 13. Comparison of the age distribution of clusters 
with masses larger than 1O 4 7 M0 and the one derived from 
the best fit model. Top: Number of clusters per age bin. 
Bottom: The cluster formation rate (number per Myr). 

been incorporated in the iV-body models that predict the 
disruption times in tidal fields are: 

1. Variations in the stellar IMF. When clusters are 
formed with a so-called top heavy IMF as is observed 
in the starburst galaxy M82 (Smith & Gallagher l2UUT)l . 
clusters will disperse much faster, since the disruption 
time depends on the number of stars in the clusters as 
idis oc N om (§EJ. Suppose the stellar IMF starts at 
1 M Q in stead of 0.1M Q , then the number of stars for 
a given cluster mass will be about a factor 10 lower. 
This will make the disruption time a factor 10 062 ~ 4 
lower. This would nicely explain the factor 5 difference 
in disruption time we observe. 

2. External perturbations. The TV-body models of 
Baumgardt & Makino ( 2003) calculated the disruption 
time of clusters in a smooth external potential from 
the host galaxy. In reality, the cluster will also expe- 
rience additional external perturbations, for example 
the encounters with molecular clouds. The clusters in 
our sample are in the inner 5 kpc of the galaxy, where 
most of the giant molecular clouds reside (Henry et al. 
2003] Kuno et al. I1995f) . The encounters with molec- 
ular clouds can speed up the disruption of clouds sig- 
nificantly (e.g. Terlevich llMTI Theuns HWT|) . 
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3. Out of equilibrium formation of clusters. All 

clusters in the TV-body models start in virial equilib- 
rium and in tidal equilibrium with the host galaxy. 
Kroupa (2004J has shown that after the gas removal 
phase, the clusters are not in virial equilibrium any- 
more and the outer parts of the cluster have expanded. 
It will be easier to dissolve these clusters than when 
all stars are in tidal equilibrium and within the tidal 
radius imposed by the host galaxy. 

4. Variations in the central concentration. The N- 
body models of Baumgardt & Makino 1|2003[1 start 
clusters with concentration values of Wo = 5-7. This 
is the average concentration of globular clusters in our 
Milky Way (Harris 1996). When clusters start with 
much smaller concentration, the core of the clusters 
is less compact and the cluster will be more vulner- 
able for external perturbations. The concentration of 
clusters in M51 can not be determined due to lack of 
resolution, but we know from young open clusters in 
the Milky Way that they have much smaller concen- 
tration indices than the globular clusters (Binney & 
Tremaine lrM7| . 

So far the clusters in M51 have not been checked for vari- 
ations in the IMF in the way that it has been done for 
clusters in other galaxies (e.g. Smith & Gallagher 120011 
Larsen et al. 120041 Maraston et al. l2004j> . Also, no A-body 
experiments have been performed including the effect of 
a tidal field and perturbations by giant molecular clouds. 
Argument 3 and 4 are based upon unknown observables 
of young clusters and they could hold for clusters in other 
galaxies as well. 

When the disruption of clusters is indeed as short as 
we derived, young massive clusters (Mj ~ 1O 6 M ) will 
not survive longer than 3.5 xlO 9 yr. This means that the 
disk of M51 is not the right location for young globular 
clusters to survive over a Hubble time. 

8. Conclusions 

We have compared the cluster population of M51 with 
theoretical predictions including evolutionary mass loss, 
cluster disruption, variable cluster formation rate and a 
magnitude limit. The age vs. mass diagrams of the ob- 
served cluster populations are binned to acquire two di- 
mensional number density plots, which can be compared 
with simulated cluster samples. The results can be sum- 
marized as follows: 

1. Artifacts introduced by our age- fitting routine do not 
systematically bias our sample towards young or old 
clusters. We present a method to correct observations 
of a cluster population for artifacts introduced by the 
age-fitting method applied. 

2. The size of the largest cluster decreases with age, from 
15 pc for clusters with log(Age/yr) < 7 to 10 pc for 
clusters with ages around 1 Gyr. In addition, the slope 
of the radius distribution seems to gets steeper in time: 



rj = 2.0 ± 0.4 for clusters younger than log(Age/yr) = 
7.5 and r\ — 2.5 ± 0.6 for clusters with log(Age/yr) > 
7.5. Both these results seem to suggest that smaller 
clusters have a larger chance to survive. However, the 
radius distribution of globular clusters in our Milky 
Way is very similar to these values found: -q = 2.4 ± 
0.2. Samples with higher spatial resolution and more 
clusters are needed to study the radius dependence. 

3. There are more old clusters at larger distances from 
the galactic center. The ratio of the number of clusters 
in the outer parts of the galaxy (3-5 kpc) over the num- 
ber of clusters in the inner part (1-3 kpc) per age bin 
increases with a factor of 1.8 in age (from log(Age/yr) 
= 6.5 to log(Age/yr) = 8.5), which is to be expected 
since the disruption time depends on the distance to 
the galactic center. 

4. Assuming that the cluster disruption time depends on 
the initial mass of the cluster as tdi s cx M? ', and us- 
ing 7 = 0.62 based on theoretical and observational 
studies (Lamers, Gieles & Portegies Zwart OHJUSJ , we 
find a typical disruption for a 10 4 M Q cluster of £4 = 

l.otoi x 108 y r > where 

we assumed a constant cluster 

formation rate. 

5. When 7 and £4 are varied together, the value found 
for 7 is similar to that predicted by Lamers, Gieles 
& Portegies Zwart (120051 based on observational and 
A-body studies. A value of 7 = 0.651°; ^5 an d ti — 
LOT 

a35 x 10 s yr are the best combination. 

6. We studied the degeneracy between formation increase 
and disruption. Models where the cluster formation 
rate increases linearly in time do not affect the disrup- 
tion timescale much (£4 gets a factor 2 higher). When 
we include bursts at the moments of encounter with 
NGC 5195, the typical disruption time is also a factor 
of 2 higher. 

7. When clusters of 1O 4 M are disrupted within 2 x 10 s 
years, and considering the power-law dependence of 
the disruption time scale with the initial cluster mass, 
even clusters with a mass of 10 6 M Q will not survive 
longer than 3.5 Gyr. This means that the disk of M51 
is not a preferred location to form a new generation of 
globular clusters. This might explain why there so far 
are no old (> Gyr) massive (> 10 6 M Q ) clusters known 
in the disks of spiral galaxies, although they are still 
forming (e.g. Westerlund 1 in the Galactic disk (Clark 
& Negueruela 2004) and the young globular cluster in 
NGC 6946 (Larsen et al. lM)T|l L 
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